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' We clarify the quantitative nature of the cosmological evolution of a global string network, that 

is, the energy density, peculiar velocity, velocity squared, Lorentz factor, formation rate of loops, 
and emission rate of Nambu-Goldstone bosons, based on a new type of numerical simulation of 
scalar fields in Eulerian meshes. We give a detailed explanation of a method to extract the above- 

| mentioned quantities to characterize string evolution by analyzing scalar fields in Eulerian meshes 

from a Lagrangian viewpoint. We confirm our previous claim that the number of long strings per 
horizon volume in a global string network is smaller than in the case of a local string network by a 

' factor of ~10 even under cosmological situations, and its reason is clarified. 

. PACS numbers: 98.80.Cq BROWN-HET-1332, OU-TAP-188 
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"sj- ■ I. INTRODUCTION 

m ; 

The idea of spontaneous symmetry breaking in high energy physics has profound implications on the vacuum 
. structure of the Universe and our Universe has presumably experienced various phase transitions, which are thermal 
or nonthermal 0- Their consequences may be traced by topological defects that may have been produced with 
them. Indeed, many implications of topological defects for cosmology have been investigated [p. Furthermore, 
recently, defect formations have also been discussed in the context of phase transitions which occur in the laboratory. 
For example, defect formations in 3 He Q and 4 He [f| are studied in detail. Thus the cosmological scenario of defect 
formation can be tested by experiments in the laboratory |6(. 

Among several types of topological defects, strings hold a unique position in cosmology, because they do not overdose 
the Universe, unlike magnetic monopoles or domain walls, settling down to a scaling solution in which the typical scale 
of the network grows in proportion to the horizon scale Q, Q • The key mechanism to achieve a scaling behavior is 
intercommutation of infinite strings to dissipate their energy by producing closed loops, which subsequently decay to 
radiation of relativistic particles or gravitational waves depending on their property. While one may understand the 
qualitative nature of the scaling solution analytically, more quantitative features such as the number of long strings 
per horizon volume or the size spectrum of loops produced cannot be obtained unless full numerical analyses are 
performed. 

Although there exist two types of strings, namely, local and global strings depending on the nature of the symmetry 
breaking, only the former have been investigated extensively for a long time as far as numerical analysis is concerned. 
By numerically solving the equation of motion of string segments derived from the Nambu-Goto action several 
groups have confirmed the scalin g be havior |sL ITfl ITTl H^. and estimated the scaling parameter £ as £ ~ 10 in the 
radiation dominated universe [HJ, [ll], C3 ■ Here, £ is defined as 

i = Poot 2 iii (i) 

where is the energy density of long strings and /i is the string tension per unit length. Thanks to this feature, 
local strings generate density fluctuations with a scale-invariant spectrum and their cosmological consequences were 
investigated extensively some time ago. Recent observations of the cosmic microwave background anisotropy [l3| . 
however, disfavor the cosmic-string scenario of structure formation, and the motivations to investigate local strings as 
a source of primordial density fluctuations have somewhat diminished, although a hybrid model of structure formation 
may still be viable, where primordial fluctuations are comprised of adiabatic fluctuations induced by inflation and 
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isocurvature perturbations by topological defects 13. 1 

On the other hand, global strings are much better motivated in the context of axion cosmology. They are formed as 
a consequence of the breaking of the Peccei-Quinn U(l) symmetry 0,0], which was introduced to solve the strong 
CP problem in quantum chromodynamics. These global strings radiate axions as associated Nambu-Goldstone (NG) 
bosons |18| . which are one of the most promising candidates for cold dark matter. Despite their importance, the 
cosmological evolution of global strings has been less studied and the results of local strings have often been borrowed 
even though there is a decisive difference between them. For local strings, the gradient energy of scalar fields is 
canceled out by gauge fields far from the core. The string core is well localized and the false vacuum energy of the 
core dominates the system. Hence the Nambu-Goto action is suitable as an effective action in order to study the 
cosmological evolution of local strings except at crossing ||| . On the other hand, for global strings, there are no gauge 
fields to cancel the gradient energy of the NG scalar field, which dominates over the false vacuum energy of the core. 
The effective action appropriate for global strings is not the Nambu-Goto action but the Kalb-Ramond action, which 
incorporates NG bosons and their couplings with the core Also, due to the gradient energy of the NG scalar 

field, a long-range force works between global strings. Thus, the behavior of the two types of string is expected to 
be different and it is nontrivial whether global strings relax to a scaling regime. In fact, for example, the behavior 
of global monopoles is quite different from that of local monopoles due to long-range forces. While the former may 
be useful in cosmology [2(J,|2l|], the latter causes a disaster [22| unless diluted by inflation. Furthermore, it has been 
shown in the literature that local and global strings behave differently in two dimensional space 0] • 

There have been several attempts to investigate the cosmological evolution of the global string network by use of 
the Kalb-Ramond action |25Ll26j . However, the Kalb-Ramond action is too complicated to be dealt with numerically. 
It has difficulty because of logarithmic divergence due to the self-energy of the string. In such a situation the authors 
and Kawasaki made the first numerical investigations of the cosmological evolution of global strings without resort 
to the Kalb-Ramond action. We instead solved the equations of motion for scalar fields forming strings in three 
dimensional Eulerian meshes j^]. We found that the global string network would also go into a scaling regime but 
the scaling parameter £ was found to be of the order of unity [27L l2^ | , which is significantly smaller than the case of 
local strings (£ ~ 10). 

Recently, however, this quantitative difference was questioned and it was claimed that the smallness of the sc aling 
parameter of global strings might be a fake due to the small dynamic range of the numerical simulations [2^, 130) . 
The authors of |29l l30l | claimed that in our simulations global strings lost their energy by excessive direct emission of 
NG bosons. Based on the speculation that such direct emission of NG bosons from long strings would be negligible 
on cosmological scales, they reached a conjecture that both types of strings behave quantitatively in the same way in 
the cosmological context, namely, £ ~ 10, with the only difference between them being the energy loss mechanism of 
closed loops. 

In order to clarify the evolution of global strings on cosmological scales, we should first study which is the dominant 
energy loss mechanism of long strings in numerical simulations, loop production or direct emission of NG bosons. We 
note that the master equation for the energy density of long strings, p^, can be expressed as 

-^p = -2H(1 + (v 2 ))poo - Tioop/Ooo - TngPoo, (2) 

where the second and the third terms on the right-hand side represent energy loss due to loop formation and direct 
emission of NG bosons or axions, respectively, and (v 2 ^ denotes the average square velocity of string segments. For 
our purpose, we need to calculate both the loop production rate ri oop and the emission rate of NG bosons Tngi an d 
compare their magnitudes in simulations. 

If the system relaxes to the scaling regime, which will be confirmed to be the case shortly, the string network is 
described by the so-called one-scale model 2 with the characteristic scale L = Poo , which grows with the horizon 
scale L oc t. If we introduce the loop production coefficient c and the emission coefficient k of the NG bosons by 

TloopPoo = C ^~~' TngPoo = (3) 



1 Note that topological defects can be compatible with cosmic inflation Q. Not e also that local strings may still be important in that 
they may emit massive particles as sources of ultrahigh energy cosmic rays ll r l . 

2 Although improvements of the sim plest on e-scale model have been proposed by taking into account effects of small-scale structures l.'ill 
and time evolution of the velocity l-'lt l32j , the original one-scale model is sufficient here because small-scale structures are smeared in 
the case of global strings and the velocity remains constant in the scaling regime, as will be shown later. 
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these parameters remain constant in the one-scale picture and are related to the scaling parameter £ as 

Indeed, if k incorrectly turned out to be much larger than c, we would find a smaller value of £ than it should really 
be. Hence, it is essential to evaluate the parameters c and k in the scaling regime. 

In our previous simulations [23] > however, it was impossible to calculate these quantities for two reasons. First, in 
the previous work a lattice was identified as part of a string based on the value of the potential energy there, which 
had the microscopic problem that we occasionally found disconnected string pieces, although the overall features 
were traced reasonably well. Second, it was impossible to monitor intercommutation of strings for lack of dynamical 
information, namely, the velocity of each string segment. These problems are overcome by our new identification 
scheme and Lagrangian analysis of the evolution of strings [33l ] . 

The rest of the paper is organized as follows. In Sec. II we present the method of our new procedure to follow 
the Lagrangian evolution of global strings in Eulerian simulations, that is, new methods of identification of strings, 
measurement of string velocity, and estimation of the intercommutation rate. In Sec. Ill the results are described 
and applied to cosmological situations. Finally, Sec. IV is devoted to a discussion and conclusion. 



II. NUMERICAL SIMULATIONS 



A. Formulation 



We consider the following Lagrangian density for two-component real scalar fields 4> a {x) (a — 1,2) which can 
accommodate global strings: 

£[<M = \g^<t>*W<t>a) - V[4> a ,T], (5) 

in the spatially flat Robertson- Walker spacetime, 

ds 2 = g^dx^dx" = dt 2 - R 2 (t)dx 2 , (6) 
with R(t) being the scale factor. We adopt the following potential at finite temperature T: 

V[cf> a ,T} = h^-^f + ^-T 2 ^ 2 (7) 
4 6 

= jW a "'/ a ) a + £(° 4 -•**)• (8) 

Here 4> 2 = <j>\ + <fe and rj 2 = a 2 - T 2 /3 = a 2 (l - T 2 /T 2 ) with T c = V3er being the critical temperature. Below the 
critical temperature, global U(l) symmetry is broken to form global strings. 
The equations of motion for the scalar fields are given by 

Mx) + 3H(t)Mx) - _L-VVa(aO + ^- = 0, (9) 

where an overdot represents the time derivative. The discretization of these differential equations is given in Appendix 
IaI Our numerical calculations are based on the staggered leapfrog method with second order accuracy both in time 
and in space. In the radiation dominated Universe, the Hubble parameter H(t) = R(t)/R{t) and cosmic time t are 
given by 

H ^ = wi^ ' = i s ^ < 10 > 

where Mpi = 1.2 x 10 19 GeV is the Plank mass and c/» is the total number of degrees of freedom for the relativistic 
particles. We define a dimensionless parameter £ as 



e / 45M 



2 \ 1/2 



C=-=(-^£L] , (11) 
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and take £ = 10 and the self-coupling A = 0.08 for definiteness, but these particular choices do not affect the results. 
We start the numerical simulation at the temperature T{ — 2T C corresponding to tj = t c /A and adopt as an initial 
condition the thermal equilibrium state with a mass equal to the inverse curvature of the potential at that time. In 
this state <f) a and </> Q are Gaussian distributed with the correlation functions 



W (^(3/)|/?> equ ^ time = / fj 1 coth ^ +m2 e^ x -^S ab , (12) 

J ( 27r ) 2Vfe +m 2 2 

, n \; i \ ; , wr,\ f d 3 k \/k 2 + m 2 (3\/k 2 +m 2 ,ikf-r-?;w ,,„\ 
(/3|0 a (a;)0 b (y)|/3) oqua i timc = J — g COt 2 ab ' ( ^ 

where m 2 = V"[4> a ,Ti] and (3 = . 4> a { x ) an d 4>a{v) are uncorrected for x ^ y. We generate the initial 
configuration in the momentum space, where the scalar fields are uncorrelated. Then they are transformed into the 
position space by fast Fourier transformation. 

We perform simulations in five different sets of lattice sizes and spacings as shown in Table Q] to investigate their 
effects on the results. In all cases, the time step is taken as At — Q.Olti, and the periodic boundary condition is 
adopted. In the case (a), the box size is nearly equal to the horizon length H" 1 and the lattice spacing to the typical 
core width S ~ 1.0/(v / 2Act) of a string at the final time tf = 200ii. The other cases have equal or larger simulation 
volumes. 

In this type of numerical calculation, it is a nontrivial task to identify string cores from the data of scalar fields 
because a point with <j) a = corresponding to a string core is not necessarily situated at a lattice point. So, in the 
next subsection, we give our new method of identification of string cores. 



B. Identification method of strings 

In previous work ^27, 28], whether a lattice point was a part of string core was judged based on the potential energy 
density there. That is, a lattice point was identified to be at the core of a string if the potential energy density there 
was found to be larger than that of a static cylindrically symmetric solution of a global string at r = Ax p h ys /V2 off 
center, where Ax p h ys is the physical length of the lattice spacing at each time. 

Although this method worked fairly well to evaluate the scaling parameter, it was inadequate to fully identify strings 
particularly for small loops which do not resemble the static cylindrically symmetric solution due to the curvature. 
As a result we occasionally found disconnected string segments that should not exist in reality. Furthermore, it was 
impossible to find a more correct position of the string core in a box beyond the lattice spacing, which is extremely 
important to evaluate the length and velocity of the string correctly. 

Here we develop a new method of string identification based on the fact that strings lie on the intersection of two 
surfaces 4>i(x, t) = and 4>2(x, t) — 0, so that if a string penetrates a (sufficiently small) plaquette, </> a has a different 
sign at one or two corners of the plaquette from the rest for each a. First, we classify the relative phase of the scalar 
fields into three groups as shown in Fig. ^ (left), that is, 

= arccos^=|L= + 27r[l-e(0 2 )], (14) 

vn + n 

« < e < |, 

. . 7T 37T 

(ii) 2~ T' (15) 

(iii) ^ < e < 2tt, 

where the arccosine should take the principal value and 0(^2) is the step function. 

Then, we judge whether a string penetrates a plaquette by monitoring the phase rotation around it just as in the 
Vachaspati-Vilenkin algorithm j3J]. If we assigned an equal range 2tt/3 of the relative phase to all regions (i)-(iii) as 
in the original Vachaspati-Vilenkin algorithm and judged the presence of a string from the phase rotation, we would 
occasionally identify a plaquette as containing a string even if <pi or c/>2 takes the same sign at its four corners. This is 
the main reason why all regions (i)-(iii) do not have equal ranges of the relative phase in our scheme, in which this 
case is avoided and 4> a takes a different sign for each a at one or two corners of each plaquette along which a circular 
phase rotation is observed. Then, using the values of ip a at its four corners we can draw two lines corresponding to 
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cj>± = and <j)2 — and their intersection is identified with the point where a string penetrates the plaquette, as shown 
in Fig.H 

In fact, the intersection point could be found outside the plaquette. In this case, the nearest point on the edge of 
the plaquette from the intersection is identified as the position of the string as shown in Fig-EJ In our simulations we 
encountered such a case very rarely and it did not cause any serious problems. 

By joining these points, the strings are completely connected and a more accurate total length of the global string 
can be obtained. One may wonder if our biased classification of the relative phase may cause some artificial effects 
on the result of numerical simulations. So we have tried another classification of the relative phase as shown in Fig.^ 
(right), that is, 

(ii) \ < < tt, (16) 

(iii) tt < 9 < 2tt, 

and confirmed that the results do not depend on the classification scheme for the relative phase. 

At each time step, we can evaluate the total length of the global string by the method shown above. Since strings 
typically move at a speed close to unity, we should multiply the length of each string segment by n — 7/i s to calculate 
the total energy density of strings. Here 7 is the Lorentz factor and [i s ~ 2nri 2 In [t/^f 1 / 2 )] 1S the line density of a 
static string |3(| . So it is important to establish a method to calculate the velocity and Lorentz factor of the string 
segments, which will be presented shortly. We can then evaluate the scaling parameter £ from Eq. 



C. Velocity of strings 

Here we describe our method to evaluate the velocity of strings, which is a nontrivial task in Eulerian calculations 
of scalar field configurations. 3 

First we expand the scalar fields <j) a {x, to + around (f> a (xo, to) up to first order, 

(j) a {x, t + St) = (f> a (x , t ) + X7cf) a (xo, t ) ■ (x - x ) + <fr a {x , t )6t (a = 1, 2). (17) 

Suppose that a string core exists at a point Xo at time to and moves to a point x at t = to + St, that is, (j) a (xo, to) = 
and 4> a (x, to + St) = for each a. Then, from Eq. H17[) we find that the loci of the string core at t = to + St lie on the 
intersection of two planes, 

V0 Q (a;o, t ) • (x ~ x ) + <fr a {x , t )St = 0, (18) 

with a — 1,2. Since motion tangential to a string is a gauge mode, we should evaluate the velocity normal to 
it. Suppose that the line normal to the string segment at (xo,to) reaches across the above intersection line at 
x = xi(xo, to, St). Then we can easily obtain the velocity of this string segment as 

/ , n xi(x ,t ,5t) - xo 

v(x ,t ) = hm 



St — »o St 

L V0 2 - 



|V<£i x V0 2 | 



(19) 

X ,t 



We calculate the velocity at each point where the strings cross plaquettes. For this purpose we need to evaluate (f> a 
and Wcf> a at arbitrary points on a plaquette, as shown explicitly in Appendix [51 where <j) a and V</> Q are given by the 
quantities on lattice points up to the second order. Collecting the value of Eq. I|19|l at each intersection of strings 
with plaquettes, we can obtain the average of the velocity, the velocity squared, and the Lorentz factor. 



3 The measurement of the velocity was also attempted in Ref. |30H . In their method, the velocity is obtained by comparing the positions 
of strings at different times. Although some efforts have been made to reduce errors, their method may cause some systematic errors 
because the time separation cannot be shorter than that sufficient to ensure a string can move to another lattice. On the other hand, 
our method is much superior; it is based on first principles and the velocity of each string segment is given by a local quantity. 
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D. Intercommutation of strings 

Now we discuss the energy loss rate from infinite strings to loops by calculating the intercommutation rate c. In 
simulations of cosmic strings based on the Nambu-Goto action, one has to assign the intercommutation probability of 
intersecting strings by hand due to the lack of microscopic information. Since the strings arc identified from the scalar 
field configuration in our simulation, we can unambiguously calculate how strings intersect with each other, namely, 
whether they simply pass through or intercommute upon collision. We can therefore identify new loops by monitoring 
the lengths of all long strings and loops at each time step and comparing the data with those at the previous time 
step. 

E. NG boson emission 

It is very difficult to directly evaluate the NG boson emission rate because separation of the emitted axions and 
NG phase associated with the string core is nontrivial. Since loops also emit these particles, it is a formidable task 
to identify axions radiated from long strings only. However, we have already described how we can obtain p^, (v 2 ), 
and ri 00 p from simulation data at each time step, and dpoo/dt can also be calculated from p^, at two adjacent time 
steps. So all the quantities in the master equation @ can be calculated from our simulation data except for the last 
term, which can be found from the equation itself. In particular, in the scaling regime, when the one-scale description 
suffices, we can find k from 

l-(v 2 ) 

k = -| 1 - c. (20) 

III. RESULTS 
A. Scaling parameter 

In Fig.01 the time evolution of £ is depicted for all cases [(a)-(e)]. The filled squares represent the time development 
of £ for the new identification method. For comparison, the results of our previous identification method based on the 
potential energy density |23,|2g, which slightly overestimates £ by a factor of 1.2, are shown by blank circles. The 
results of the original Vachaspati-Vilenkin identification scheme are indicated by blank squares, obtained by counting 
the number of boxes through which a string passes as was done in Ref . [3(j • This method overestimates £ by a factor 
of 1.4. 

We find that differences in the simulation settings except the box size normalized by the final horizon scale do not 
affect the overall behavior significantly. £ becomes smaller in the smallest box size with = H^ 1 at the final time. This 
is mainly because under the periodic boundary conditions, a string feels an attractive force from the boundary and 
tends to disappear. Indeed, in the case (a) this artificial effect starts to operate before the system really relaxes to 
the scaling solution, so £ continues decreasing without a plateau. On the other hand, in the case (b), the boundary 
effect becomes operative after t ~ 130 when £ starts to decrease from the scaling value. The other three cases, whose 
box size is larger than 2if _1 , are free from the boundary effects and the scaling parameter £ remains constant. Hence 
we conclude that after the relaxation period the global string network goes into the scaling regime. From Table ITU 
which shows average values of £ after some relaxation period (t > 80), £ is found to be £ ~ 0.80 in the scaling regime. 

B. String velocity 

In Figs. E]-[7| the time evolutions of the average velocity (v), the average square velocity (v 2 ), and the Lorentz 
factor (7) are shown. After some relaxation period, they become fairly constant in all cases. The average values are 
given in Table[H] We find that the average values become larger in the smallest box size. This is mainly because under 
the periodic boundary conditions, a string feels an attractive force from the boundary and tends to be accelerated. 
However, such an effect becomes negligible for a box size larger than as observed in Table [H] Thus we find 

the average velocity (v) ~ 0.60 and the average square velocity (v 2 ) ~ 0.50 » (v) 2 , in the scaling regime. On the 
other hand, the average Lorentz factor has a large scatter in time, although the long-time average is fairly constant 
with (7) ~ 1.8. Since string segments moving with a speed close to the light velocity have extremely large Lorentz 
factors and push up the average dramatically, the fluctuation in the number of such string segments results in this 
large scatter. This is the reason the average Lorentz factor (7) = 1.8 ± 0.2 is much larger than 1/yl — (v 2 ) and 
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l/Jl - (v). Thus the energy per unit length of a string is enhanced by a factor (7) ~ 1.8 compared with a static 
string. 

C. Energy dissipation coefficients 

As explained in Sec. Ill Dl we can obtain the intercommutation rate c without assigning its probability by hand. 
The result is given in Table |H] where we find that c is larger in the smallest box size. This can be understood in the 
same way, namely, under the periodic boundary conditions, a string feels an attractive force from the boundary and 
tends to intercommute more often. Since such an effect is unimportant for a box size larger than 2H~ 1 as discussed 
above, c is found to be c = 0.40 ± 0.04. 

Taking account of the relation (|20() , the emission parameter n is given in Table [H] In the simulations with the 
smallest box, strings intercommute more often due to the boundary effect, which suppresses NG boson emission. But 
again such an effect is inoperative for cases with box size larger than 2H~ 1 . So we conclude that the NG boson 
emission rate is n = 0.16 ± 0.04. 



IV. DISCUSSION AND CONCLUSION 

In this paper we have investigated the cosmological evolution of the global string network in detail. In our nu- 
merical simulations, the equations of motion of the two-component real scalar field are solved on Eulerian meshes. 
In order to follow the time evolution of global strings, we have developed a new identification method which enables 
us to find a more correct position of the string core in a box beyond the lattice spacing. Furthermore, we have 
given a detailed explanation of the method to extract Lagrangian quantities, such as velocity and intercommutation, 
characterizing the evolution of global strings. The NG boson emission rate is obtained from the master equation 
Thus the quantitative nature of the cosmological evolution of the global string network is elucidated without setting 
the intercommutation probability of two intersecting strings by hand. Specifically, we find the scaling parameter 
characterizing the energy density £ ~ 0.80, the peculiar velocity (v) ~ 0.60, the velocity squared (v 2 ) ~ 0.50, the 
Lorentz factor (7) = 1.8 ± 0.2, the formation rate of loops c = 0.40 ± 0.04, and the emission rate of Nambu-Goldstone 
bosons k — 0.16 ± 0.04 from our results. 

In our simulations, due to the limit of the dynamic range of the lattices, the logarithm of the ratio of the Hubble 
radius to the string width, which appears in the expression for the effective line density of global strings, took 
ln(t/<$) ~ 5, while in the cosmological setting it can be as large as 0(100). The authors of |29l l30| argue that k is 
inversely proportional to \xi(t/S), although c and (i> 2 ) do not have such dependence. Based on this speculation, they 
claimed that we obtained a smaller value of the scaling parameter £ in our previous work [2]} than it should really 
be, because k was incorrectly large there. They also argued that in the cosmological situation loop production is the 
dominant mechanism of energy dissipation of long strings and that £ would take the same value as for local strings. 

Now that we have all the values of relevant quantities from our simulation data, we can disprove their claim. Indeed, 
we find that k is smaller than c and hence direct emission of NG bosons does not dominate the energy loss mechanism 
even in our setting of numerical simulations with a relatively small dynamic range, ]xi(t/S). If we extrapolate the value 
of k with the scaling k oc 1/ ln(t/<5), we find n < 0.01 in the cosmological situation. Even then the scaling parameter 
(0} does not increase much, yielding £ = 1.6 ± 0.3, which is still much smaller than that of local strings. 

Thus we conclude that there is a quantitative difference between the cosmological evolution of global strings and that 
of local strings. It is based on the qualitative difference that global strings have a long-range force and intercommute 
more often with larger (i> 2 ) than do local strings. Note that it is by no means surprising that global defects behave 
differently from local defects in cosmology. In the case of monopoles, as discussed in the Introduction, global monopoles 
evolve in a strikingly different manner than do magnetic monopoles. In the case of strings, both local and global 
defects relax to a scaling solution and their difference is more subtle. So it is not until our numerical analysis of the 
evolution of the latter from the Lagrangian point of view is performed that their difference is fully elucidated. 

Finally, we use our results to obtain a constraint on the symmetry breaking scale 77 = f a of the Peccei-Quinn U(l) 
symmetry. From £ = 1.6 ± 0.3 and (7) = 1.8 ± 0.2, we find f a < (0.16 - 1.2) x 10 12 GeV for the normalized Hubble 
parameter h = 0.7 27]. We also note that the Lagrangian method developed in this paper is directly applicable to 
other species of extended objects, for instance, topological defects such as local strings and nontopological solitons 
such as Q balls as well. 
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Note added in proof 

After we submitted the original manuscript, we became aware of an analytic estimate of £ which reports that it is 
~ 1 even in the case of the local string network. We are grateful to M. Hindmarsh for informing us of his result |35| . 
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APPENDIX A: DISCRETIZATION OF DIFFERENTIAL EQUATIONS 

The equation of motion of the scalar fields is given by 

fax) + 3H(t)Mx) - J^V 2 M*) + = 0. (Al) 
In order to discretize the above equation, it is reduced to first-order differential equations, 



fa — "a 



1 „o , , s dV 



7T a = -3H(t)7T a -j^V 2 Mx) 



Expanding <p a (t, x) and <f> a {t + At, x) around the intermediate time step t + \At, we find 



(A2) 



<j> a {t + At,x)-4> a {t,x) = Atj> a (t+^At,x\ + 0((Atf), (A3) 

up to the third order. Thus 7r a at the intermediate time step is given by 

n a (t + \ At, „) = + - ^ X) + 0{{Atf) . (A4) 

In the same way, 7T a and n a at the time step t are represented by the quantities at the two adjacent intermediate steps 
t±\At, 

= " aft + ^' 3;) At 7rQft "^ a;) + 0((A^), 

«•„(*, x) = ^+k^)+^-h^, x ) +0{{Mn (A5) 

The second-order derivatives are approximated up to the second order in A as 

-0(A 2 ), 
0{A 2 ), 

0{A 2 ). (A6) 





i>a{t,x) 


(f> a (t,x + A,y,z) 


- 2(f) a (t,x,y, z) + 


4> a (t,x- A,y,z) 




dx 2 




A 2 




d 2 . 


i>a{t, x ) 


</>a(t,x,y + A,z) 


- 2<f> a (t,x,y, z) + 


(j) a (t,x,y- A, z) 




dy 2 




A 2 




d 2 . 




cj> a (t,x,y,z + A) 


- 2cj) a (t,x,y, z) + 


cj) a (t,x,y,z - A) 




dz 2 




A 2 





Thus the fundamental equations are discretized up to the second order in both space and time, 
Ux) + M{t)Ux) - — L_vVaW + dV g^ X)] 
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b a (t + iAt, x) - 4> a {t - ^At, x) 4> a {t + §At, sc) + a (t - iAt, 03 
\-oH{t) 



R(t) 2 ^ 



At w 2 

■ a (t,x + At)-2</> a (t,x)+<f> a (t,x- Ai)\ , 0V[0„(t,a:)] 



A 2 



30a 



0((Ai) 2 ,A 2 



= 0, 



(A7) 



where A,, = (A, 0,0), A y = (0, A,0), and A z = (0,0, A). 
In summary, in our numerical simulations 



t+ i^At.x 



1 



1 + 3H(t)At/2 
+At 



3H(t)At\ 



>„ [t--At,x 



— y 



2 J 

b a {t, x + A t )- 2cf> a (t, x) + cf> a (t, x - AO 1 dV[<t> a (t, x)} 



<j) a {t + At,x) = <t> a {t,x) + At(j) a ( t+ -At, a; 



(A8) 



The value of a at the time step f, which is required to calculate the string velocity, is evaluated from Eq. (|A5 



APPENDIX B: QUANTITIES AT AN ARBITRARY POINT ON A PLAQUETTE 

In order to evaluate the velocity of a string correctly, we need to evaluate quantities at an arbitrary point on a 
plaquette. That is, <j) a and V0 a within a plaquette should be expressed by their values at its four corners. As an 
example, let us consider a plaquette parallel to the z plane with four corners (x, y, z), [x + A, y, z), (x, y + A, z), and 
(x + A, y + A, z) and express <j> a {t, x + a A, y + j3 A, z) and V^ (t, x + a A, y + j3A, z), which we denote collectively by 
II(x + aA, y + (3 A) below, in terms of their values at these four points. Here < a < 1 and < (3 < 1. 

We express II at the four corners by an expansion around (x + aA, y + (3 A) as 

U.(x,y) = n{x + aA,y + pA)-aA^-pA^ 

+ r A V + a/3A + ^ A V + ° (A ); 

an an 

n{x + A,y) = IL(x + aA,y+j3A) + (l-a)A— - PA- 
ox ay 

- °» 2a2 - <* - + ?' ,2&2 S + °< a3 >- 

n(x,y + A) = n(.T + aA,y + /3A)-aA^ + (l-/3)A^ 

ox oy 

n{x + A,y + A) = n{x + aA,y + /3A) + {l-a)A^- + (l-P)A^- 

ox oy 

+5(1 - «)^ 2 + P - "K 1 - « A2 S + 5 (1 - « 2A2 + ° (A3) - (B1) 

Making an appropriate combination, Tl(x + aA, y + j3A) can be expressed by its values at the four vertices in the 
plaquette, 



(1 - a)(l - (3) n{x, y) + a(l - /3) H(i + A, y) + (1 - a)/9 E(a, y + A) + a/3 n(ai + A, y + A) 

= n(x + aA, y + PA) + ~a(l - «)A 2 + |/3(1 - /?)A 2 + G(A 3 ) (B2) 
= n(a; + aA,y + /3A)+0(A 2 ). (B3) 
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In particular, replacing II(x + aA, y + (3 A) by <j> a (x + aA, y + /3A), it can be expressed by the quantities on the lattice 
points up to second order, 



(j> a (t, x + aA,y + (3A,z) = (1 - a)(l - 0) Q (i, x, y, z) + a(l - 0) j> a (t,x + A,y, z) 

+ (1 - a)(3 o (i, x, y + A, z) + a(3 4> a (t, x + A, y + A, z) + C(A 2 ). 

In the same way, inserting d<f> a / dx into II, we find 

<f> a (t, x + aA, y + (3 A, z) 



dx 



ft ut a\ d ^S-> x ^V^ z ) i /, a \ 9<l>a(t,x + A,y,z) 



+(l-a)/3 



<9x 

9^ a (t, x,y + A,z) 



a/3 ■ 



dx 

b a (t,x + A,y + A,z) 



1 

2A 



dx dx 
(1 - o)(l - 13) {0 Q (t, x + A, y, z) - <j> a {t, x- A,y, z)| 

+a(l - /3) |?!> a (t,a;+ 2A, y, z) - a (i, x, y, z) j 

+ (1 - a)/3 {<p a (t, x + A, y + A, z) - (f> a (t, x — A,y + A, z) j 

+a/3 {^ a (t, x + 2A, y + A, z) - a (t, x, y + A, z)} 



0(A 2 ). 



Here we have used the relations 

(/) a (x + A, y, z)-(j) a (x- A,y, z) _ dcf) a (x,y,z) 



2A 

(j) a {x + 2A, y, z) - (j) a (x, y, z) 
2A 

Interchanging x for y and a for (3, d<p a / dy is also given by 



dx ■ °< a2 > 

a (x + A,y,z) 



dx 



0{A Z ). 



cj) a (t,x + aA,y + (3A,z) dcf> a (t,x,y,z) dcf> a (t,x + A,y,z) 
= (1 — ct)(l— p) q \- a(l — p) — 



dy 



dy 



dy 



1 

2A 



dz 



dy dy 
(1 - a)(l - (3) \j> a {t, x, y + A, z) - o (t, x, y - A, z)} 

- /?) x + A, y + A, z) - a (t, x + A, y - A, z)} 

+ (1 - a)(3 x, y + 2A, z) - a (i, x, y, z)| 

+a/3 |^ B (t, x + A, y + 2A, z) - <0 B (t, x + A, y, z)} + C(A 2 ). 



0(A 2 ) 



Using the relations obtained above and Eq. (|B2(1 . d(f> a /dz is calculated as 
</> a 0,x + oA,y + /3A,z) _ ^ a (t, x + aA, y + /3A, z + A) - <f> a (t, x + aA, y + (3A, z - A) 



2A 



1 

2A 



(1 - a)(l - /?) x, y, z + A) - a (i, x, y, z - A)} 

+a(l - /?) |0 o (f, x + A, y, z + A) - 4> a (t, x + A, y, z - A)} 
+ (1 - a)(3 {^ a (t, x, y + A, z + A) - a (i, x, y + A, z - A)} 

+a/9 {0 B (i, x + A, y + A, z + A) - o (t, x + A, y + A, z - A)} 



(B4) 



G(A 2 ) (B5) 



(B6) 



(B7) 



_ dMt,x,V + \z) + ap ^ a (t,x + A,y + A,z) + c(a9) (Bg) 



(B9) 



-— a(l — a 
4 



v , ( d 2 4> a (t, x + aA, y + /3A, z + A) _ <9 2 a ft, x + aA, y + (3 A, z - A) \ 
' \ dx 2 dx 2 J 
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7/3(1 ~/3)A 



d 2 <j) a (t, x + a A, y + f3A, z + A) d 2 cj) a {t 1 x + a A, y + /3A, 2 



<9y 2 



4" 

+0(A 2 ). 

Although these approximations appear of first-order accuracy, using the relations, 
d 2 <j> a (t, x + aA, y + f3A, z + A) d 2 cj) a (t, x + aA, y + (3 A, z - A) 



dy 2 



dx 2 



2A 



9a; 2 

d 3 4> a (t,x + aA,y + f3A,z) 
dx 2 dz 



+ 0(A 3 ), 



d 2 cj) a (t, x + aA, y + f3A, z + A) d 2 </> a (*, x + aA, y + (3 A, z - A) 



dy 2 



2A 



dy 2 

d 3 (j) a {t,x + aA,y + f3A,z) 
dy 2 dz 



G(A 3 ), 



we find they have second-order accuracy. As a result, up to the second order, we find 
4> a {t, x + aA, y + (3 A, z) 



dz 



1 

2A 



(1 - a)(l - 0) {0 a (t,x, y, z + A) - a (i,x, y,z - A) I 
+a(l - 13) x + A, y, z + A) - <f> a (t, x + A, y, z - A)} 

+ (1 - a)l3 x, y + A, z + A) - <j> a (t, x, y + A, z - A)} 

+a/3 |0 o (t, x + A, y + A, z + A) - </> a (t, x + A, y + A, z - A)} 



-G(A 2 ). 
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TABLE I: Five different sets of simulations. 



Case Lattice Lattice spacing (A) £ Realization Box size/ H 





number 


[units of UR(t)] 






(at final time) 


(a) 


128 3 


V3/8 


10 


100 


1 (at 200) 


(b) 


256 3 




10 


10 


1 (at 200) 


(c) 


128 3 


\/3/4 


10 


10 


2 (at 200) 


(d) 


256 3 


V3/S 


10 


10 


2 (at 200) 


(e) 


256 3 


73/4 


10 


10 


4 (at 200) 



TABLE II: Results of numerical simulations. 

Case ^ (v) (f 2 ) 7 c k 

(a) 0.72 0.63 0.50 2.0 ± 0.2 0.52 ± 0.05 0.08 ±0.05 

(b) 0.79 0.65 0.50 1.9 ±0.3 0.48 ± 0.04 0.08 ± 0.04 

(c) 0.77 0.60 0.50 1.8 ±0.2 0.40 ± 0.02 0.17 ±0.02 

(d) 0.80 0.60 0.49 1.8 ±0.2 0.42 ± 0.03 0.16 ±0.03 

(e) 0.80 0.60 0.50 1.8 ±0.2 0.40 ± 0.04 0.16 ±0.04 




FIG. 1: Left: Relative phase of the scalar fields is classified into three groups, (i) < 6 < f , (ii) f < 6 < (iii) ^ < 6 < 2tt. 
Right: Another classification of the relative phase, (i) < 6 < ^, (ii) -| < 6 < n, (iii) tt < 9 < 2n, was also tried but the 
numerical results did not depend on these choices. 



14 





\ 
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\ 

\ 

\ 



q 2 

FIG. 2: p a and q a are points corresponding to <f> a = for each a. These points are obtained by linear interpolation using the 
values of <f> a at two corners of the plaquette between which it changes sign. The line <f> a = for each a is drawn by simply 
connecting p a and q a by a straight line. The intersection of these two lines is identified as the position through which a string 
penetrates a plaquette. 




FIG. 3: When the intersection is found outside the plaquette, the nearest point on the edge of the plaquette is identified as the 
penetration point of the string. 
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FIG. 4: Time evolution of the scaling parameter £. Filled squares represent our new identification method. Blank circles 
correspond to our previous identification method based on the potential energy density [27I l2Sjl . Blank squares are results of 
the Vachaspati-Vilenkin algorithm. 
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FIG. 5: Time evolution of average velocity of global strings is shown. 
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FIG. 6: Time evolution of average square velocity. 
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FIG. 7: Time evolution of average Lorentz factor. 



